!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>

!>\brief
!! This module provides subroutines for the definition of global
!! degrees of freedom for continuous finite element spaces.
!!
!! \details
!!
!! The definition of continuous finite element spaces is somewhat more
!! complex that the definition of discontinuous spaces, due to the
!! continuity constraint. In fact, while in the construction of
!! discontinuous finite element spaces it suffices to consider
!! vertexes (\f$0\f$-simplexes), sides (\f$(d-1)\f$-simplexes) and
!! elements (\f$d\f$-simplexes), as discussed in \c mod_grid, when
!! constructing continuous spaces one must ensure continuity on each
!! family of \f$e\f$-simplexes for \f$e=0,\ldots,d\f$. For instance,
!! for three dimensional problems, one also needs to consider the
!! element edges (\f$1\f$-simplexes), which are never used in three
!! dimensional discontinuous elements.
!!
!! In this module, the nodes of the continuous finite element space
!! are placed and numbered as follows:
!! <ul>
!!  <li> each vertex corresponds to one node, and these nodes are
!!  numbered as the corresponding vertices
!!  <li> nodes defined on the reference \f$e\f$-simplex are mapped on
!!  each \f$e\f$-dimensional face of the grid for \f$e=1,\ldots,d\f$,
!!  and these nodes are numbered in increasing order following the
!!  numbering of the \f$e\f$-dimensional faces.
!! </ul>
!! Notice that \f$(d-1)\f$-faces are in fact the sides of the grid,
!! but the numbering used for the definition of such faces in this
!! module can be different from the numbering used in \c mod_grid.
!!
!! \section ddc_dofs Domain decomposition
!!
!! The domain decomposition information is handled in a similar way to
!! what is done for the grid vertexes. In particular, to each degree
!! of freedom is attached information concerning:
!! <ul>
!!  <li> its global "identity"
!!  <li> the neighbouring subdomains sharing it.
!! </ul>
!!
!! \subsection ghost_dofs Dofs associated with ghost entities
!!
!! In \c mod_bcs, some ghost entities are defined. The degrees of
!! freedom associated with such entities are defined as ghost dofs,
!! since they can be useful for finite element formulations with a
!! nonstandard stencil, such as the side based stabilizations. The
!! ghost dofs carry a limited amount of information and are attached
!! to the ghost entities. As local dofs, ghost dofs are unique, while
!! as global dofs they can be repeated. More precisely, global dofs
!! which belongs to more than one ghost entity appear multiple times
!! as ghost dofs, each time with a different local index; however, all
!! the local copies share the same global index.
!!
!! \subsection ddc_dofs_MPI MPI usage
!!
!! The MPI buffers are handled as follows (all the deallocations are
!! done in \c new_cgdofs):
!! <table>
!!  <tr><td></td> <th>allocate</th><th>set</th></tr>
!!  <tr><td>ds_send_buf  </td>
!!            <td>new_cgdofs      </td> <td>insert_face     </td></tr>
!!  <tr><td>ds_recv_buf  </td>
!!            <td>new_cgdofs      </td> <td>receive_ddc_data</td></tr>
!!  <tr><td>idx_send_buf </td>
!!            <td>new_cgdofs      </td> <td>clear_face      </td></tr>
!!  <tr><td>send_buf     </td>
!!            <td>new_cgdofs      </td> <td>clear_face      </td></tr>
!!  <tr><td>recv_buf     </td>
!!            <td>receive_ddc_data</td> <td>receive_ddc_data</td></tr>
!!  <tr><td>send_conf_buf</td>
!!            <td>insert_ddc_data </td> <td>insert_ddc_data </td></tr>
!!  <tr><td>recv_conf_buf</td>
!!            <td>send_ddc_data   </td> <td>send_ddc_data   </td></tr>
!! </table>
!! \note This module uses various nonblocking MPI communications.
!! According to the discussion found <a
!! href="https://svn.mpi-forum.org/trac/mpi-forum-web/ticket/143">here</a>
!! and <a
!! href="https://svn.mpi-forum.org/trac/mpi-forum-web/ticket/47">here</a>,
!! we avoid any use of <tt>mpi_request_free</tt>, and terminate each
!! send or receive operation with a wait operation.
!<----------------------------------------------------------------------
module mod_cgdofs

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_mpi_utils, only: &
   mod_mpi_utils_initialized, &
   mpi_integer, wp_mpi,   &
   mpi_status_size,       &
   mpi_irecv, mpi_isend,  &
   mpi_request_null,      &
   mpi_wait, mpi_waitall, &
   mpi_bcast

 use mod_octave_io, only: &
   mod_octave_io_initialized, &
   write_octave, &
   read_octave,  &
   read_octave_al

 use mod_perms, only: &
   mod_perms_initialized, &
   t_perm, t_dperm, &
   operator(.eq.),  &
   operator(.lt.),  &
   operator(**),    &
   perm_table,      &
   idx,             &
   dperm_reduce

 use mod_master_el, only: &
   mod_master_el_initialized, &
   t_me, t_lagnodes, &
   new_me, clear

 use mod_grid, only: &
   mod_grid_initialized, &
   t_ddc_gv, t_ddc_nv,  &
   t_grid, t_ddc_grid,  &
   affmap, write_octave

 use mod_bcs, only: &
   mod_bcs_initialized, &
   t_bcs, b_dir

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_cgdofs_constructor, &
   mod_cgdofs_destructor,  &
   mod_cgdofs_initialized, &
   t_cgdofs, t_ddc_cgdofs, &
   t_ghost_dof,          &
   new_cgdofs,           &
   clear,                &
   write_octave

 private

!-----------------------------------------------------------------------

! Module types and parameters

 !> Collection of global degrees of freedom.
 !!
 !! The fields \c dofs and \c x are ordered as follows:
 !! <ul>
 !!  <li> the first <tt>ndofsd(0)</tt> columns coincides with the grid
 !!  vertices;
 !!  <li> the following <tt>ndofs(1)</tt> columns are placed on
 !!  \f$1\f$-dimensional faces
 !!  <li> \f$\ldots\f$
 !!  <li> the last <tt>ndofs(d)</tt> columns are internal degrees of
 !!  freedom.
 !! </ul>
 !! Notice that this ordering can be exploited to identify which nodes
 !! can be eliminated in the static condensation.
 type t_cgdofs
   integer :: d, m
   !> total number of dofs
   integer :: ndofs
   !> total number of dofs on \f$e\f$-dimensional faces, ndofsd(0:d)
   integer, allocatable :: ndofsd(:)
   !> the \c dofs table maps local degrees of freedom of element \c ie
   !! to global degrees of freedom as <tt>dofs(:,ie)</tt>
   integer, allocatable :: dofs(:,:)
   real(wp), allocatable :: x(:,:) !< node coordinates (m,ndofs)
   !> Dirichlet boundary degrees of freedom.
   !!
   !! This field has one column for each Dirichlet degree of freedom,
   !! and three rows indicating:
   !! <ol>
   !!  <li> the index of the Dirichlet boundary degree of freedom
   !!  <li> the value of the type parameter \c btype (see \c mod_bcs)
   !!  <li> the value of the boundary region \c breg (see \c mod_bcs)
   !! </ol>
   !!
   !! A degree of freedom is of Dirichlet type if it belongs to an
   !! \f$e\f$-face of the grid which vertices are all of Dirichlet
   !! type; the type parameter and the boundary region index are
   !! chosen as the maximum value among all the vertices of the face.
   !! These choices are consistent with the way boundary conditions
   !! are transferred from boundary sides to boundary vertices in \c
   !! mod_bcs.
   !! 
   !! \c dir_dofs is (partially) ordered according to \c btype, and
   !! the degrees of freedom on the boundary region \c ib are found in
   !! <tt>dir_dofs(:,idir_dofs(1,ib):idir_dofs(2,ib))</tt>.
   integer, allocatable :: dir_dofs(:,:)
   integer, allocatable :: idir_dofs(:,:) !< used to index \c dir_dofs.
   !> Natural degrees of freedom
   !!
   !! This field collects the indices of all the degrees of freedom
   !! which are not of Dirichlet type, i.e. all the degrees of freedom
   !! which do not appear in \c dir_dofs. This includes internal
   !! degrees of freedom and boundary degrees of freedom with a
   !! Neumann boundary condition. The last <tt>ndofsd(d)</tt> entries
   !! of <tt>nat_dofs</tt> correspond to the nodes collocated in the
   !! element interior: these are the degrees of freedom eliminated in
   !! a static condensation process.  The two arrays of indexes
   !! <tt>dir_dofs(1,:)</tt> and \c nat_dofs allow a very
   !! straightforward implementation of the essential boundary
   !! conditions.
   integer, allocatable :: nat_dofs(:)
 end type t_cgdofs

 !> Ghost degree of freedom
 !!
 !! The local index of a ghost dof is
 !! <ul>
 !!  <li> consistent with the local dof numbering if the dof already
 !!  belongs to the local subdomain
 !!  <li> larger than that of any local dof if the dof does not belong
 !!  to the local subdomain.
 !! </ul>
 type t_ghost_dof
   integer :: i          !< local index
   integer :: i_dir_nat  !< local index as Dirichlet or natural dof
   integer :: gi         !< global index
   logical :: dir        !< whether the dof is of Dirichlet type
   integer :: gi_dir_nat !< global index as Dirichlet or natural dof
   real(wp), allocatable :: x(:) !< coordinates
   !> the following fields are significant only for dir dofs
   integer :: btype
   integer :: breg
 end type t_ghost_dof

 !> Domain decomposition information.
 !!
 !! The ddc information stored for each dof is analogous to the ddc
 !! information carried by vertices. In particular, it is of two
 !! types: global indexing and indexing among neighbouring subdomains.
 !!
 !! \note The global ordering of the degrees of freedom has similar
 !! properties to the local ordering. More precisely:
 !! <ul>
 !!  <li> the first <tt>ddc_grid\%ngv</tt> degrees of freedom coincide
 !!  with the vertexes
 !!  <li> the following degrees of freedom are collocated on the
 !!  element faces (without any further specification)
 !!  <li> the last <tt>lagnodes(grid\%d)\%nn*dc_grid\%nge</tt> degrees
 !!  of freedom are internal to the elements.
 !! </ul>
 !! This is very important, because it ensures that the global indexes
 !! of vertex and face degrees of freedom do not change after a static
 !! condensation of the element degrees of freedom.
 !!
 !! \note Concerning \c gnat_dofs, we ensure that the internal dofs
 !! are numbered after the vertex and face dofs, again to allow a
 !! simple implementation of the static condensation.
 type t_ddc_cgdofs
   !> total number of global dofs, same value on all subdomains
   integer :: ngdofs
   !> total number of global dofs on \f$e\f$-dimensional faces, same
   !! value on all subdomains, ndofsd(0:d)
   integer, allocatable :: ngdofsd(:)
   integer :: nndofs !< # neighbour dofs
   !> number of neighbour dofs for each domain (0:nd-1)
   integer, allocatable :: nndofs_id(:)
   type(t_ddc_gv), allocatable :: gdofs(:) !< global dofs (ndofs)
   type(t_ddc_nv), allocatable :: ndofs(:) !< neighbour dofs (nndofs)
   integer :: ngdir_dofs !< global number of dir dofs
   integer :: ngnat_dofs !< global number of nat dofs
   !> global numbering of the Dirichlet degrees of freedom
   integer, allocatable :: gdir_dofs(:)
   !> global numbering of the natural degrees of freedom
   integer, allocatable :: gnat_dofs(:)
   !> ghost element dofs (size(dofs,1),nns)
   type(t_ghost_dof), allocatable :: ghost_dofs(:,:)
   !> number of ghost dofs (i.e. number of dofs which belong to ghost
   !! entities and <em>do not belong to the local subdomain</em>:
   !! local dofs connected also to ghost entities are not included)
   integer :: n_ghost_dofs
   integer :: n_ghost_dir_dofs !< analogous to \c n_ghost_dofs
   integer :: n_ghost_nat_dofs !< analogous to \c n_ghost_dofs
   !> indexes in \c ghost_dofs of the ghost dofs which do not belong
   !! to the local subdomain
   integer, allocatable :: added_ghost_dofs(:,:)
 end type t_ddc_cgdofs

 !> Type for the side linked list
 type t_face_list
   type(t_dperm) :: v               !< vertices
   integer, allocatable  :: in(:)   !< nodes
   real(wp), allocatable :: xn(:,:) !< node coordinates
   logical :: dir = .false. !< face of Dirichlet type
   integer :: btype         !< Dirichlet boundary type
   integer :: breg          !< Dirichlet region
   type(t_face_list), pointer :: next => null()
   !> The following fields are specific for ddc
   logical :: ddc = .false.
   integer, allocatable :: ddc_cand(:,:)  !< candidate subdomains
   type(t_ddc_nv), allocatable :: ndofs(:)!< neighbour dofs
   integer, allocatable :: ing(:)         !< global dofs indexes
 end type t_face_list

! Module variables

 ! public members
 logical, protected ::               &
   mod_cgdofs_initialized = .false.
 ! private members
 character(len=*), parameter :: &
   this_mod_name = 'mod_cgdofs'

 ! MPI buffers
 integer :: &
   ids, ide, & !< start and end indexes: 0, ddc_grid%nd-1
   ils, ile, & !< start and end indexes for the lower \id buffers
   ihs, ihe    !< start and end indexes for the higher \id buffers
 integer, target, allocatable :: &
   !> <tt>ds_send_buf(ddc_grid%id+1:ddc_grid%nd-1)</tt>: for each
   !! subdomain, indicates the number of columns in \c idx_send_buf,
   !! send_buf and \c recv_conf_buf (and the same for
   !! <tt>ds_recv_buf</tt>)
   !!
   !! \c ds_send_buf: allocation
   ds_send_buf(:), ds_send_req(:), &
   ds_recv_buf(:), ds_recv_req(:), &
   !> <tt>idx_send_buf(:,idx)</tt> has the following structure:
   !! <ul>
   !!  <li> first element: domain index
   !!  <li> second element: face dimensiion \f$e\f$
   !!  <li> following \c nn elements: dof indexes
   !! </ul>
   idx_send_buf(:,:), &
   !> <tt>send_buf(:,idx)</tt> and <tt>recv_buf(:,:)</tt> have the
   !! following structure:
   !! <ul>
   !!  <li> first element: the dimension \c e of the face
   !!  <li> following \f$e+1\f$ elements: face vertexes, with the
   !!  indexing of the receiver and the ordering of the sender
   !!  <li> following \c nn elements: dof indexes of the sender
   !!  <li> following \c nn elements: global dof indexes (decided by
   !!  the sender)
   !! </ul>
   send_buf(:,:), send_req(:), offset(:), &
   recv_buf(:,:), &
   !> <tt>send_conf_buf(:,idx)</tt> has the following structure:
   !! <ul>
   !!  <li> first element: 1 if used, zero otherwise
   !!  <li> following \c nn elements: the local indexes on the other
   !!  subdomain
   !! </ul>
   send_conf_buf(:,:), send_conf_req(:), &
   recv_conf_buf(:,:), recv_conf_req(:)

 interface clear
   module procedure clear_cgdofs, clear_ddc_cgdofs
 end interface

 interface write_octave
   module procedure write_cgdofs, write_ddc_cgdofs, &
     write_ghost_dof_struct
 end interface write_octave

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_cgdofs_constructor()
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.)  .or. &
       (mod_kinds_initialized.eqv..false.)     .or. &
       (mod_mpi_utils_initialized.eqv..false.) .or. &
       (mod_octave_io_initialized.eqv..false.) .or. &
       (mod_perms_initialized.eqv..false.)     .or. &
       (mod_master_el_initialized.eqv..false.) .or. &
       (mod_grid_initialized.eqv..false.)      .or. &
       (mod_bcs_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_cgdofs_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif

   mod_cgdofs_initialized = .true.
 end subroutine mod_cgdofs_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_cgdofs_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_cgdofs_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_cgdofs_initialized = .false.
 end subroutine mod_cgdofs_destructor

!-----------------------------------------------------------------------

 !> Definition of the continuous finite element degrees of freedom
 !! (i.e. Lagrangian nodes). Since most of the work is represented by
 !! the definition of the grid \f$e\f$-faces for \f$e=0,\ldots,d\f$,
 !! whenever one needs many finite element spaces on the same grid
 !! it's convenient build them all together, hence the use of vector
 !! arguments.
 !! \bug The arguments of \c new_cgdofs should be vectors.
 !!
 !! The algorithm used is very similar to what is done in \c mod_grid
 !! for the construction of the side list.
 !!
 !! Concerning the domain decomposition information, we proceed as
 !! follows:
 !! <ul>
 !!  <li> each time a new face is created, we check whether it can be
 !!  shared and find out which are the subdomains which are
 !!  <em>candidates</em> for sharing it by looking at the intersection
 !!  of the sharing subdomains for all the face vertexes (the key
 !!  point here is that a subdomain can share all the vertexes of a
 !!  face without sharing the face itself, hence the need to find
 !!  candidates and then verify a posteriori their feedback)
 !!  <li> after defining the face list, the data sent by the
 !!  subdomains with lower \c id are received (this is done in two
 !!  steps: first receive the size of the data and then the data
 !!  themselves), included in the face list and then a confirmation
 !!  message is sent back to inform about the amount of used data
 !!  (feedback)
 !!  <li> then the local degrees of freedom are stored back in the \c
 !!  t_cgdofs object and the face list is dismantled, and while doing
 !!  this the send buffers are prepared to send the data to the
 !!  subdomains with larger \c id
 !!  <li> finally, the data are sent to the subdomains with larger \c
 !!  id, and the feedback is received and used to complete the \c
 !!  t_ddc_cgdofs object.
 !! </ul>
 !! \note Each subdomain sends its information to all the candidate
 !! subdomains, however, since these latter don't know that they
 !! are being seen as candidates from another subdomains, the
 !! communication must be done in two steps:
 !! <ol>
 !!  <li> send the size of the data
 !!  <li> send the data.
 !! </ol>
 !<
 subroutine new_cgdofs(dofs,lagnodes,grid,bcs , comm,ddc_grid,ddc_dofs)
  type(t_lagnodes), intent(in)  :: lagnodes(:)
  type(t_grid),     intent(in)  :: grid
  type(t_bcs),      intent(in)  :: bcs
  type(t_cgdofs),   intent(out) :: dofs
  !> optional ddc specific arguments
  integer,            intent(in),  optional :: comm
  type(t_ddc_grid),   intent(in),  optional :: ddc_grid
  type(t_ddc_cgdofs), intent(out), optional :: ddc_dofs

  integer :: e, iv, ie, ifl, i_shift1, i_shift2
  integer :: ndofsd_el(0:grid%d)
  integer, allocatable :: ii(:), nat_dofs(:)
  type(t_dperm) :: face_list_tail
  type(t_face_list), allocatable, target :: face_list(:,:)
  type(t_me) :: medm1
  character(len=*), parameter :: this_sub_name = 'new_cgdofs'
  ! ddc specifi variables
  logical :: ddc
  !> Counter for the creation of the global degrees of freedom. The
  !! vertex degrees of freedom are numbered as the corresponding the
  !! vertexes, for the remaining dofs, each subdomain adds its own
  !! dofs counted by this variable. Counting the degrees of freedom
  !! separately according to the face dimension allows us to place the
  !! global element dofs after the global vertex and face dofs.
  integer, allocatable :: gdof_dcounter(:)
  !> This variable is used as a temporary for ddc_dofs%ndofs, since
  !! the correct size for such array is known only at the end of the
  !! subroutine.
  type(t_ddc_nv), allocatable :: n_dofs(:)
  integer :: ierr
  character(len=1000) :: msg
   
   !--------------------------------------------------------------------
   ! Set here some fields which are easily computed
   dofs%d = grid%d
   dofs%m = grid%m
   dofs%ndofs = grid%nv ! will be used and changed later
   allocate(dofs%ndofsd(0:grid%d));
   dofs%ndofsd(0) = grid%nv
   dofs%ndofsd(1:grid%d-1) = 0 ! will be changed later
   dofs%ndofsd(  grid%d  ) = lagnodes(grid%d)%nn * grid%ne
   ! number of dofs for each element
   ndofsd_el(0) = grid%d+1                    ! vertices
   do e=1,grid%d-1                            ! faces
     ndofsd_el(e) = lagnodes(e)%nn * grid%me%children(e)%ns
   enddo
   ndofsd_el(grid%d) = lagnodes(grid%d)%nn    ! elements
   allocate(dofs%dofs(sum(ndofsd_el),grid%ne))
   allocate(dofs%idir_dofs(2,bcs%nbreg))

   ddc = present(ddc_grid)
   if(ddc) then
     ids = 0;             ide = ddc_grid%nd-1
     ils = 0;             ile = ddc_grid%id-1
     ihs = ddc_grid%id+1; ihe = ddc_grid%nd-1
     allocate( ds_send_buf(ihs:ihe) , ds_send_req(ihs:ihe) )
     allocate( ds_recv_buf(ils:ile) , ds_recv_req(ils:ile) )
     ds_send_buf = 0; ds_recv_buf = 0
     ! The remaining MPI buffers can not be allocated now because
     ! their size must first be determined.
     allocate(gdof_dcounter(0:grid%d))
     gdof_dcounter(0)      = ddc_grid%ngv
     ! Notice that
     ! a) the number of faces is not know a priori
     ! b) gdof_dcounter(grid%d) will be used as a counter for the
     !   creation of the global element dofs
   else
     ! These indexes are initialized since they are used in the
     ! definition of some argument bounds
     ids = 0; ide = 0
     ils = 0; ile = 0
     ihs = 0; ihe = 0
   endif
   !--------------------------------------------------------------------

   !--------------------------------------------------------------------
   ! The total number of nodes is not easily determined beforehand, so
   ! it's better to start from the face nodes, which are those
   ! requiring the largest part of the work, and then adding nodes
   ! placed on vertices and internal nodes, which definition is
   ! trivial.
   dofs%idir_dofs(1,:) = 1; dofs%idir_dofs(2,:) = 0
   ! allocate the face list
   allocate(face_list(grid%nv,grid%d-1))

   do e=1,grid%d-1
     if(lagnodes(e)%nn.gt.0) then

       ! terminate the lists with grid%nv+1, larger than any iv
       face_list_tail = dperm_reduce((/(grid%nv+1, ifl=1,e+1)/))
       do iv=1,grid%nv
         allocate(face_list(iv,e)%next)
         face_list(iv,e)%next%v = face_list_tail
       enddo

       ! build the list and set dofs%dofs
       do ie=1,grid%ne ! element loop
         do ifl=1,grid%me%children(e)%ns ! local face loop
           call insert_face(grid,face_list(:,e), dofs%ndofs,        &
                       dofs%ndofsd(e),dofs%dofs(:,ie),              &
                       grid%e(ie)%iv(grid%me%children(e)%s(:,ifl)), &
                       e , ifl , lagnodes(e) , ndofsd_el , ddc_grid,&
                       ds_send_buf)
         enddo
       enddo

     endif
   enddo
   ! complete the list with the Dirichlet bcs
   call new_me(medm1,grid%d-1)
   do iv=1,bcs%nb_side ! boundary side loop
     if(bcs%b_side(iv)%bc.eq.b_dir) then ! set to dir the faces
       do e=1,grid%d-2
         if(lagnodes(e)%nn.gt.0) then
           do ifl=1,medm1%children(e)%ns ! local face loop
             call setdir_face1(face_list(:,e),                      &
                   bcs%b_side(iv)%s%iv(medm1%children(e)%s(:,ifl)), &
                   bcs%b_side(iv)%btype, -bcs%b_side(iv)%s%ie(2)    )
           enddo
         endif
       enddo
       e = grid%d-1 ! the side itself, not included in medm1%children
       if(lagnodes(e)%nn.gt.0) then
         call setdir_face1(face_list(:,e),                          &
                   bcs%b_side(iv)%s%iv, & ! use all the side vertexes
                   bcs%b_side(iv)%btype, -bcs%b_side(iv)%s%ie(2)    )
       endif
     endif
   enddo
   call clear(medm1)
   call setdir_face2( dofs%idir_dofs,grid,lagnodes,face_list,bcs%nbreg )

   ! get the ddc information if required
   if(ddc) then
     call receive_ddc_data( grid,ddc_grid,comm,lagnodes,gdof_dcounter )
     call insert_ddc_data( ddc_grid,lagnodes,face_list,ierr,msg, &
                ds_recv_buf,recv_buf,send_conf_buf,send_conf_req )
     if(ierr.ne.0) call error(this_sub_name,this_mod_name,msg)
     call send_feedback(comm)
   endif
   !--------------------------------------------------------------------

   !--------------------------------------------------------------------
   ! Add the vertex and internal nodes; clean the list retrieving also
   ! the information concerning dofs%x
   dofs%ndofs = dofs%ndofs + dofs%ndofsd(grid%d)
   allocate( dofs%x(grid%m,dofs%ndofs) )
   do iv=1,bcs%nb_vert ! count the Dirichlet vertices
     if(bcs%b_vert(iv)%bc.eq.b_dir) then
       dofs%idir_dofs(1,bcs%b_vert(iv)%breg+1:bcs%nbreg) = &
         dofs%idir_dofs(1,bcs%b_vert(iv)%breg+1:bcs%nbreg) + 1
       dofs%idir_dofs(2,bcs%b_vert(iv)%breg:bcs%nbreg) = &
         dofs%idir_dofs(2,bcs%b_vert(iv)%breg:bcs%nbreg) + 1
     endif
   enddo
   allocate( dofs%dir_dofs(3, dofs%idir_dofs(2,bcs%nbreg))       , &
             dofs%nat_dofs(dofs%ndofs-dofs%idir_dofs(2,bcs%nbreg)) )
   if(ddc) then
     allocate( ddc_dofs%gdofs(dofs%ndofs) )
     ddc_dofs%nndofs = 0 ! will be incremented
     allocate( ddc_dofs%nndofs_id(ids:ide) ); 
     ddc_dofs%nndofs_id = 0 ! will be incremented
     allocate( n_dofs(dofs%ndofs) ) ! work array
     ! we biuld also some working arrays
     iv = 2+maxval(lagnodes%nn) ! used as temporary
     ie = sum(ds_send_buf)
     allocate(idx_send_buf(iv,ie))
     iv = 1+maxval( (/(e+1 , e=1,grid%d-1)/) + 2*lagnodes%nn )
     allocate(send_buf(iv,ie)); send_buf = -1
     allocate(send_req(ihs:ihe))
   endif
   
   ! vertices
   do ie=1,grid%ne
     dofs%dofs(1:grid%d+1,ie) = grid%e(ie)%iv
   enddo
   do iv=1,grid%nv
     dofs%x(:,iv) = grid%v(iv)%x
   enddo
   ! here the second row of idir_dofs is reset and used as counter
   dofs%idir_dofs(2,:) = dofs%idir_dofs(1,:)-1
   do iv=1,bcs%nb_vert
     if(bcs%b_vert(iv)%bc.eq.b_dir) then
       ifl = dofs%idir_dofs(2,bcs%b_vert(iv)%breg) + 1 ! as temporary
       dofs%dir_dofs(1,ifl) = bcs%b_vert(iv)%v%i
       dofs%dir_dofs(2,ifl) = bcs%b_vert(iv)%btype
       dofs%dir_dofs(3,ifl) = bcs%b_vert(iv)%breg
       dofs%idir_dofs(2,bcs%b_vert(iv)%breg) = ifl
     endif
   enddo
   if(ddc) then
     ! vertexes already have a global indexing: the simplest thing
     ! here it to number the vertex degrees of freedom after the
     ! vertexes
     do iv=1,grid%nv
       ddc_dofs%gdofs(iv)%gi = ddc_grid%gv(iv)%gi
       if(ddc_grid%gv(iv)%ni.gt.0) then ! shared dof
         ddc_dofs%nndofs = ddc_dofs%nndofs+1
         ifl = ddc_grid%gv(iv)%ni ! used as temporary
         ddc_dofs%nndofs_id(ddc_grid%nv(ifl)%id) = &
           ddc_dofs%nndofs_id(ddc_grid%nv(ifl)%id) + 1
         ddc_dofs%gdofs(iv)%ni = ddc_dofs%nndofs
         ! this is simply a copy of the neighbour vertex (this happens
         ! because vertex dofs are the first ones both locally and
         ! globally)
         n_dofs(ddc_dofs%nndofs) = ddc_grid%nv(ifl)
       endif
     enddo
   endif

   ! faces (and clean the lists)
   if(ddc) then
     allocate( offset(ihs:ihe) )
     offset = 0 ! correct also when offset is empty
     do e=ihs+1,ihe
       offset(e) = offset(e-1) + ds_send_buf(e-1)
     enddo
   endif
   do e=1,grid%d-1
     if(lagnodes(e)%nn.gt.0) then
       do iv=1,grid%nv
         ! in case of domain decomposition, this also fills
         ! idx_send_buf and send_buf
         call clear_face( face_list(iv,e) , dofs%x , dofs%idir_dofs ,&
                  dofs%dir_dofs , e , lagnodes(e)%nn , ddc_grid ,    &
                  gdof_dcounter , grid%d , ddc_dofs , n_dofs ,       &
                  idx_send_buf , send_buf , offset )
       enddo
     endif
   enddo
   if(ddc) deallocate(offset)
   deallocate(face_list)

   ! elements
   if(lagnodes(grid%d)%nn.gt.0) then
     allocate(ii(lagnodes(grid%d)%nn))
     ii = (/(ifl, ifl=1,lagnodes(grid%d)%nn)/)
     i_shift1 = sum(dofs%ndofsd(0:grid%d-1)) ! non-internal nodes
     do ie=1,grid%ne
       i_shift2 = i_shift1 + (ie-1)*lagnodes(grid%d)%nn
       dofs%dofs(sum(ndofsd_el(0:grid%d-1))+ii,ie) = i_shift2 + ii
       dofs%x(:,i_shift2+ii) = affmap(grid%e(ie),lagnodes(grid%d)%x)
     enddo
     deallocate(ii)
     if(ddc) then
       do iv=i_shift1+1,dofs%ndofs
         ! Here we number the element dofs without considering the
         ! vertex and face dofs. The numbering will be corrected
         ! afterwords, when we know the total number of global face
         ! dofs.
         gdof_dcounter(grid%d) = gdof_dcounter(grid%d) + 1
         ddc_dofs%gdofs(iv)%gi = gdof_dcounter(grid%d)
       enddo
     endif
   endif

   ! fix the natural degrees of freedom
   allocate(nat_dofs(dofs%ndofs))
   nat_dofs = (/ (ifl, ifl=1,dofs%ndofs) /)
   nat_dofs(dofs%dir_dofs(1,:)) = -1
   dofs%nat_dofs = pack(nat_dofs,nat_dofs.gt.0)
   deallocate(nat_dofs)
   !--------------------------------------------------------------------

   !--------------------------------------------------------------------
   ! Complete the ddc data using the feedback of the neighboring
   ! subdomains.
   if(ddc) then

     ! send data to the candidates and receive the feedback
     call send_ddc_data(grid,ddc_grid,comm,lagnodes,gdof_dcounter)
     ddc_dofs%ngdofs  = sum(gdof_dcounter)
     allocate(ddc_dofs%ngdofsd(0:grid%d))
     ddc_dofs%ngdofsd = gdof_dcounter
     ! Fix the global numbering of element dofs
     i_shift1 = sum(dofs%ndofsd(0:grid%d-1)) ! non-internal nodes
     i_shift2 = sum(gdof_dcounter(0:grid%d-1))
     ddc_dofs%gdofs(i_shift1+1:dofs%ndofs)%gi =            &
         ddc_dofs%gdofs(i_shift1+1:dofs%ndofs)%gi + i_shift2

     ! complete the neighbour information
     call insert_feedback_data(ddc_dofs,n_dofs,lagnodes, &
                               idx_send_buf,recv_conf_buf)

     deallocate( ds_send_buf, ds_send_req, ds_recv_buf, ds_recv_req, &
                 idx_send_buf, send_buf, send_req, recv_buf,         &
                 send_conf_buf, send_conf_req, recv_conf_buf,        &
                 recv_conf_req )

     ! trim the neighbour dofs array
     allocate(ddc_dofs%ndofs(ddc_dofs%nndofs))
     ddc_dofs%ndofs = n_dofs(1:ddc_dofs%nndofs)
     deallocate(n_dofs)

     ! define the global numbering of dir and nat dofs
     call set_g_dir_nat_data(ddc_dofs,dofs,ddc_grid,comm)

     ! define the ghost dofs
     call set_ghost_dofs(ddc_dofs,grid,dofs,ddc_grid,comm)

   endif
   !--------------------------------------------------------------------

 end subroutine new_cgdofs

!-----------------------------------------------------------------------

 pure subroutine insert_face( grid , face_list , ndofs , ndofsd , dofs,&
     iv , e , ifl , lagnodes , ndofsd_el , ddc_grid , ds_sb )
  type(t_grid),      intent(in) :: grid
  integer,           intent(in) :: iv(:), e, ifl, ndofsd_el(0:)
  type(t_lagnodes),  intent(in) :: lagnodes
  type(t_face_list), intent(inout), target :: face_list(:)
  integer,           intent(inout) :: ndofs, ndofsd, dofs(:)
  type(t_ddc_grid),  intent(in), optional :: ddc_grid
  integer,           intent(inout) :: ds_sb(ihs:)

  integer :: ii(lagnodes%nn), i_shift, i
  real(wp) :: b(grid%m,e)
  type(t_dperm) :: iv_red
  type(t_face_list), pointer :: p, new_face
     
   iv_red = dperm_reduce(iv) ! reduce to basic form
   ii = (/ (i, i=1,lagnodes%nn) /)
   i_shift = sum(ndofsd_el(0:e-1)) + (ifl-1)*lagnodes%nn
   ! insert in the linked list
   p => face_list(iv_red%x(1)) ! the first vertex idexes the array
   search: do
     if(iv_red.lt.p%next%v) then ! add new face to the list
       ! build the new object
       allocate(new_face)
       new_face%v = iv_red
       allocate(new_face%in(lagnodes%nn))
       new_face%in = ndofs + ii
       allocate(new_face%xn(grid%m,lagnodes%nn))
       do i=1,e
         b(:,i) = grid%v(iv_red%x(i+1))%x - grid%v(iv_red%x(1))%x
       enddo
       new_face%xn = matmul(b,lagnodes%x) &
                    + spread(grid%v(iv_red%x(1))%x,2,lagnodes%nn)
       ! Dirichlet bcs will be considered elsewhere
       ndofs  = ndofs  + lagnodes%nn
       ndofsd = ndofsd + lagnodes%nn
       ! update dofs information, use the right order!!
       dofs(i_shift+ii) = &
         new_face%in( lagnodes%stab( idx( iv_red%pi**(-1) ) , : ) )

       ! check whether the face must be sent to neighbour domains
       if(present(ddc_grid)) then
         call set_candidates( new_face%ddc_cand,ddc_grid,iv_red%x )
         ds_sb(new_face%ddc_cand(1,:)) =   &
           ds_sb(new_face%ddc_cand(1,:)) + 1
       endif

       ! insert the new object in the list
       new_face%next => p%next
       p%next        => new_face

       exit search
     elseif(iv_red.eq.p%next%v) then ! found an existing face
       ! update dofs information, use the right order!!
       dofs(i_shift+ii) = &
         p%next%in( lagnodes%stab( idx( iv_red%pi**(-1) ) , : ) )
       exit search
     endif

     p => p%next
   enddo search

 contains

   pure subroutine set_candidates(can,ddc_grid,iv)
    integer, allocatable, intent(out) :: can(:,:)
    type(t_ddc_grid), intent(in) :: ddc_grid
    integer, intent(in) :: iv(:)

    logical :: shared
    integer :: n, i, id, j, k
    integer, allocatable :: tmp(:,:)
    type(t_ddc_nv) :: nv(size(iv))

     if(any(ddc_grid%gv(iv)%ni.eq.0)) then 
       ! at least one vertex not shared
       allocate(can(1+size(iv),0))
     else
       nv = ddc_grid%nv(ddc_grid%gv(iv)%ni)
       allocate(tmp(1+size(nv),nv(1)%nd))
       n = 1
       do i=1,nv(1)%nd
         id = nv(1)%id(i) ! neighbour
         if(id.gt.ddc_grid%id) then ! candidates have larger id
           tmp(1,n) = id; tmp(2,n) = nv(1)%in(i)
           check_do: do j=2,size(nv)
             shared = .false.
             locate_do: do k=1,nv(j)%nd
               if(nv(j)%id(k).eq.id) then
                 tmp(1+j,n) = nv(j)%in(k)
                 shared = .true.
                 exit locate_do
               endif
             enddo locate_do
             if(.not.shared) exit check_do
           enddo check_do
           if(shared) n = n+1
         endif
       enddo
       n = n-1
       allocate(can(1+size(nv),n))
       can = tmp(:,1:n)
       deallocate(tmp)
     endif
   end subroutine set_candidates

 end subroutine insert_face

!-----------------------------------------------------------------------

 !> Set the dir information of the face dofs, step 1
 pure subroutine setdir_face1( face_list , iv , btype , breg )
  type(t_face_list), intent(inout), target :: face_list(:)
  integer, intent(in) :: iv(:), btype, breg

  type(t_dperm) :: iv_red
  type(t_face_list), pointer :: p

   iv_red = dperm_reduce(iv) ! reduce to basic form
   ! find the face in the linked list (we know that it exists!)
   p => face_list(iv_red%x(1)) ! the first vertex idexes the array
   search: do
     if(iv_red.eq.p%next%v) then ! found
       if(p%next%dir) then ! already known as dir face
         p%next%btype = max( p%next%btype , btype )
         p%next%breg  = max( p%next%breg  , breg  )
       else
         p%next%dir = .true.
         p%next%btype = btype
         p%next%breg  = breg
       endif
       exit search
     endif

     p => p%next
   enddo search

 end subroutine setdir_face1

!-----------------------------------------------------------------------

 !> Set the dir information of the face dofs, step 2
 !!
 !! This step requires that each face has already a well defined breg:
 !! this is the reason why it is separated from \c setdir_face1.
 pure subroutine setdir_face2( idir_dofs , grid , lagnodes , &
                               face_list , nbreg )
  type(t_grid),      intent(in) :: grid
  type(t_lagnodes),  intent(in) :: lagnodes(:)
  !> The intent is inout to be able to use pointers; otherwise it
  !! should be an intent(in) argument.
  type(t_face_list), intent(inout), target :: face_list(:,:)
  integer,           intent(in) :: nbreg
  integer,           intent(inout) :: idir_dofs(:,:)

  integer :: e, iv
  type(t_face_list), pointer :: p

   ! We loop on the whole face_list and count the dir dofs
   do e=1,grid%d-1
     if(lagnodes(e)%nn.gt.0) then
       do iv=1,grid%nv
         p => face_list(iv,e)%next
         do
          if(.not.associated(p%next)) exit ! p points to the terminator
           if(p%dir) then ! dir face
             idir_dofs(1,p%breg+1:nbreg) = &
               idir_dofs(1,p%breg+1:nbreg) + lagnodes(e)%nn
             idir_dofs(2,p%breg  :nbreg) = &
               idir_dofs(2,p%breg  :nbreg) + lagnodes(e)%nn
           endif
           p => p%next
         enddo
       enddo
     endif
   enddo

 end subroutine setdir_face2

!-----------------------------------------------------------------------

 !> While inserting the received data, we also prepare the feedback
 !! buffers
 !<
 pure subroutine insert_ddc_data( ddc_grid , lagnodes , face_list , &
                         ierr , msg , ds_rb , rb, conf_sb , conf_sr )
  type(t_ddc_grid),   intent(in)            :: ddc_grid
  type(t_lagnodes),   intent(in)            :: lagnodes(:)
  type(t_face_list),  intent(inout), target :: face_list(:,:)
  integer,            intent(out)   :: ierr
  character(len=*),   intent(out)   :: msg
  integer,            intent(in)    :: ds_rb(ils:), rb(:,:)
  integer,allocatable,intent(inout) :: conf_sb(:,:), conf_sr(:)

  integer :: dim1, id, inn_id, inn, e, i_perm, i,            &
    i_tmp_n(maxval(lagnodes%nn)), i_tmp_g(maxval(lagnodes%nn))
  type(t_dperm) :: iv_r
  type(t_face_list), pointer :: p
  character(len=26) :: formt

   dim1 = 1+maxval(lagnodes%nn)
   allocate( conf_sb(dim1,size(rb,2)) )
   conf_sb(1,:) = 0 ! initialize to not used
   allocate( conf_sr(ils:ile) )

   ierr = 0
   main_do: do id=ils,ile; do inn_id=1,ds_rb(id)
     inn = sum(ds_rb(ils:id-1)) + inn_id
     e = rb(1,inn)
     iv_r = dperm_reduce(rb(2:e+2,inn))
     p => face_list(iv_r%x(1),e)%next
     search: do
       ! It is possible that the received face does not exist: this
       ! happens when the candidate subdomain (i.e. this subdomain)
       ! does not share the face.
       if(.not.associated(p%next)) exit search ! p points to the terminator
       if(iv_r.eq.p%v) then ! found
         ! reorder the received data
         i_perm = idx(iv_r%pi)
         i_tmp_n(1:lagnodes(e)%nn) =                           &
           rb(e+2+               lagnodes(e)%stab(i_perm,:),inn)
         i_tmp_g(1:lagnodes(e)%nn) =                           &
           rb(e+2+lagnodes(e)%nn+lagnodes(e)%stab(i_perm,:),inn)
         ! update the face list
         if(.not.p%ddc) then ! new
           p%ddc = .true.
           allocate(p%ndofs(size(p%in)))
           allocate(p%ing  (size(p%in)))
           do i=1,size(p%in)
             p%ndofs(i)%i = p%in(i)
             p%ndofs(i)%nd = 1
             allocate(p%ndofs(i)%id(1)); p%ndofs(i)%id(1) = id
             allocate(p%ndofs(i)%in(1)); p%ndofs(i)%in(1) = i_tmp_n(i)
             p%ing(i) = i_tmp_g(i)
           enddo
         else ! already known neighbour dof
           ! consistency check
           if(any(p%ing.ne.i_tmp_g(1:lagnodes(e)%nn))) then
             ierr = 1
             write(formt,'(a,i3,a,i3,a)') &
             '(a,i5,a,i5,a,',size(p%ing),'i8,a,',lagnodes(e)%nn,'i8)'
             write(msg,formt) &
               'Inconsistency on proc ', ddc_grid%id, &
               ' reading data from proc ', id,        &
               '; conflicting global indexes ',p%ing, &
               ' and ',i_tmp_g
             exit main_do
           endif
           do i=1,size(p%in)
             p%ndofs(i)%nd = p%ndofs(i)%nd + 1
             call add_alloc(p%ndofs(i)%id,id)
             call add_alloc(p%ndofs(i)%in,i_tmp_n(i))
           enddo
         endif
         ! update the conf_sb
         i_perm = idx(iv_r%pi**(-1))
         conf_sb(1,inn) = 1
         conf_sb(2:1+lagnodes(e)%nn,inn) =                &
                           p%in(lagnodes(e)%stab(i_perm,:))
        exit search
       endif
       p => p%next
     enddo search
   enddo; enddo main_do

 end subroutine insert_ddc_data

!-----------------------------------------------------------------------

 pure subroutine clear_face( face_list,dofsx,idir_dofs,dofsdir_dofs, &
         e,nn,ddc_grid,gdof_dcounter,d,ddc_dofs,n_dofs , idx_sb,sb,ofs )
  type(t_face_list), intent(inout), target :: face_list
  real(wp),          intent(inout)         :: dofsx(:,:)
  integer,           intent(inout) :: idir_dofs(:,:), dofsdir_dofs(:,:)
  ! ddc specific arguments
  integer,            intent(in)              :: e, nn
  type(t_ddc_grid),   intent(in), optional    :: ddc_grid
  integer,            intent(inout)           :: gdof_dcounter(0:)
  integer,            intent(in)              :: d
  type(t_ddc_cgdofs), intent(inout), optional :: ddc_dofs
  type(t_ddc_nv),     intent(inout) :: n_dofs(:)
  integer,            intent(out)   :: idx_sb(:,:), sb(:,:)
  integer,            intent(inout) :: ofs(ihs:)

  integer :: i
  integer, allocatable :: ii(:)
  type(t_face_list), pointer :: p, p2
  ! ddc specific variables
  logical :: ddc_check
  integer :: id, j, is, ie

   ddc_check = present(ddc_grid)

   p => face_list%next
   do
    if(.not.associated(p%next)) exit ! p points to the terminator
     dofsx(:,p%in) = p%xn
     ddc_if: if(ddc_check) then
       if(p%ddc) then ! global dof already defined
         do i=1,size(p%in)
           ddc_dofs%gdofs(p%in(i))%gi = p%ing(i)
           ddc_dofs%nndofs = ddc_dofs%nndofs + 1
           ddc_dofs%nndofs_id(p%ndofs(i)%id) = &
              ddc_dofs%nndofs_id(p%ndofs(i)%id) + 1
           ddc_dofs%gdofs(p%in(i))%ni = ddc_dofs%nndofs
           n_dofs(ddc_dofs%nndofs) = p%ndofs(i)
         enddo
       else ! define a new global dof (so far not shared)
         do i=1,size(p%in)
           gdof_dcounter(e) = gdof_dcounter(e) + 1
           ddc_dofs%gdofs(p%in(i))%gi = sum(gdof_dcounter(0:d-1))
         enddo
       endif
       ! prepare the send buffer for the candidates
       do i=1,size(p%ddc_cand,2)
         id = p%ddc_cand(1,i) ! candidate index
         ofs(id) = ofs(id) + 1; j = ofs(id)
         idx_sb(1    ,j) = id
         idx_sb(2    ,j) = e
         is = 3;    ie = is-1 + size(p%in)
         idx_sb(is:ie,j) = p%in
         sb(1    ,j) = e
         is = 2;    ie = is-1 + e+1
         sb(is:ie,j) = p%ddc_cand(2:e+2,i)
         is = ie+1; ie = is-1 + nn
         sb(is:ie,j) = p%in
         is = ie+1; ie = is-1 + nn
         sb(is:ie,j) = ddc_dofs%gdofs(p%in)%gi 
       enddo
     endif ddc_if
     if(p%dir) then
       allocate(ii(size(p%in)))
       ii = idir_dofs(2,p%breg) + (/(i, i=1,size(p%in))/)
       dofsdir_dofs(1,ii) = p%in
       dofsdir_dofs(2,ii) = p%btype
       dofsdir_dofs(3,ii) = p%breg
       idir_dofs(2,p%breg) = idir_dofs(2,p%breg) + size(p%in)
       deallocate(ii)
     endif
     p2 => p
     p => p%next
     deallocate(p2)
   enddo
   ! deallocate the terminator of the list
   deallocate(p)
 end subroutine clear_face

!-----------------------------------------------------------------------

 pure subroutine insert_feedback_data(ddc_dofs,n_dofs,lagnodes, &
                                      idx_sb,conf_rb)
  type(t_lagnodes),   intent(in) :: lagnodes(:)
  integer,            intent(in) :: idx_sb(:,:), conf_rb(:,:)
  type(t_ddc_cgdofs), intent(inout) :: ddc_dofs
  type(t_ddc_nv),     intent(inout) :: n_dofs(:)

  integer :: inn, id, e, i, idof, nndofs

   do inn=1,size(conf_rb,2)
     if(conf_rb(1,inn).eq.1) then ! face was used
       id = idx_sb(1,inn)
       e  = idx_sb(2,inn)
       do i=1,lagnodes(e)%nn
         idof = idx_sb(2+i,inn)
         if(ddc_dofs%gdofs(idof)%ni.eq.0) then ! new
           ddc_dofs%nndofs        = ddc_dofs%nndofs        + 1
           ddc_dofs%nndofs_id(id) = ddc_dofs%nndofs_id(id) + 1
           nndofs = ddc_dofs%nndofs
           ddc_dofs%gdofs(idof)%ni= nndofs
           n_dofs(nndofs)%i = idof
           n_dofs(nndofs)%nd = 1
           allocate(n_dofs(nndofs)%id(1)); n_dofs(nndofs)%id(1) = id
           allocate(n_dofs(nndofs)%in(1)); n_dofs(nndofs)%in(1) =  &
                                              conf_rb(1+i,inn)
         else ! already known neighbour dof
           ddc_dofs%nndofs_id(id) = ddc_dofs%nndofs_id(id) + 1
           nndofs = ddc_dofs%gdofs(idof)%ni
           n_dofs(nndofs)%nd = n_dofs(nndofs)%nd + 1
           call add_alloc(n_dofs(nndofs)%id,id)
           call add_alloc(n_dofs(nndofs)%in,conf_rb(1+i,inn))
         endif
       enddo
     endif
   enddo

 end subroutine insert_feedback_data

!-----------------------------------------------------------------------

 !> Receive ddc data from lower \c id subdomains
 subroutine receive_ddc_data( grid,ddc_grid,comm,lagnodes, &
                              gdof_dcounter   )
  type(t_grid),     intent(in)  :: grid
  type(t_ddc_grid), intent(in)  :: ddc_grid
  integer,          intent(in)  :: comm
  type(t_lagnodes), intent(in)  :: lagnodes(:)
  ! Additional output is collected in ds_recv_buf and recv_buf
  integer,          intent(inout) :: gdof_dcounter(0:)

  integer :: i, ierr, dim1, e, is, gdof(1:grid%d), req
  integer, allocatable :: recv_sta(:,:)

   !1) Get the amount of received data
   do i=ils,ile
     call mpi_irecv( ds_recv_buf(i) , 1 , mpi_integer , &
                     i, 1, comm, ds_recv_req(i), ierr   )
   enddo
   allocate( recv_sta(mpi_status_size,ils:ile) )
   call mpi_waitall(size(ds_recv_req),ds_recv_req,recv_sta,ierr)

   !2) Receive the data
   dim1 = 1+maxval( (/(e+1 , e=1,grid%d-1)/) + 2*lagnodes%nn )
   allocate(recv_buf(dim1,sum(ds_recv_buf))) ! correct also for id.eq.0

   ds_recv_req = mpi_request_null
   is = 1
   do i=ils,ile
     if(ds_recv_buf(i).gt.0) &
       call mpi_irecv( recv_buf(1,is) , dim1*ds_recv_buf(i) ,       &
                       mpi_integer, i, 2, comm, ds_recv_req(i), ierr)
     if(i.eq.ddc_grid%id-1) call mpi_irecv( gdof(1) , grid%d , &
                       mpi_integer, i, 3, comm, req, ierr )
     is = is + ds_recv_buf(i)
   enddo
   call mpi_waitall(size(ds_recv_req),ds_recv_req,recv_sta,ierr)

   ! On all processors, gdof_dcounter(0) has been set already
   if(ddc_grid%id.ne.0) then
     call mpi_wait(req,recv_sta(:,0),ierr)
     gdof_dcounter(1:grid%d) = gdof
   else
     gdof_dcounter(1:grid%d) = 0
   endif

   deallocate(recv_sta)

 end subroutine receive_ddc_data

!-----------------------------------------------------------------------

 subroutine send_feedback(comm)
  integer, intent(in) :: comm

  integer :: dim1, is, i, ierr

   dim1 = size(send_conf_buf,1)
   is = 1
   send_conf_req = mpi_request_null
   do i=ils,ile
     if(ds_recv_buf(i).gt.0) &
       call mpi_isend( send_conf_buf(1,is) , dim1*ds_recv_buf(i),     &
                       mpi_integer, i, 4, comm, send_conf_req(i), ierr)
     is = is + ds_recv_buf(i)
   enddo
   ! The wait are at the end of the whole cgdofs creation,
   ! because we don't need these data any more.
   
 end subroutine send_feedback

!-----------------------------------------------------------------------

 !> Send ddc data to lower \c id subdomains
 !!
 !! Send the ddc data and start the receives of the confirmation
 !! messages.
 !<
 subroutine send_ddc_data(grid,ddc_grid,comm,lagnodes,gdof_dcounter)
  type(t_grid),     intent(in)    :: grid
  type(t_ddc_grid), intent(in)    :: ddc_grid
  integer,          intent(in)    :: comm
  type(t_lagnodes), intent(in)    :: lagnodes(:)
  integer,          intent(inout) :: gdof_dcounter(0:)

  integer :: dim1_1, dim1_2, is, i, ierr, gdof(1:grid%d), req
  integer, allocatable :: recv_sta(:,:)

   dim1_1 = size(send_buf,1)
   dim1_2 = 1+maxval(lagnodes%nn)
   allocate( recv_conf_buf(dim1_2,size(send_buf,2)) )
   allocate( recv_conf_req(ihs:ihe) )
   gdof = gdof_dcounter(1:grid%d)

   send_req = mpi_request_null
   recv_conf_req = mpi_request_null
   is = 1
   do i=ihs,ihe

     call mpi_isend( ds_send_buf(i) , 1 , mpi_integer , &
                     i, 1, comm, ds_send_req(i), ierr )

     if(ds_send_buf(i).gt.0) &
       call mpi_isend( send_buf(1,is) , dim1_1*ds_send_buf(i) ,   &
                       mpi_integer, i, 2, comm, send_req(i), ierr )

     if(i.eq.ddc_grid%id+1) call mpi_isend( gdof(1) , grid%d , &
                       mpi_integer, i, 3, comm, req, ierr )

     ! Start receiving the feedback
     if(ds_send_buf(i).gt.0) &
       call mpi_irecv( recv_conf_buf(1,is) , dim1_2*ds_send_buf(i) ,   &
                       mpi_integer, i, 4, comm, recv_conf_req(i), ierr )

     is = is + ds_send_buf(i)
   enddo
   
   ! receive the feedback
   allocate( recv_sta(mpi_status_size,ddc_grid%id+1:ddc_grid%nd-1) )
   call mpi_waitall(size(recv_conf_req),recv_conf_req,recv_sta,ierr)

   ! Logically, this could be omitted, since once we receive the
   ! feedback we know that the send has been completed, however MPI
   ! must know that the send requests can be freed.
   call mpi_waitall(size(ds_send_req  ),ds_send_req  ,recv_sta,ierr)
   call mpi_waitall(size(send_req     ),send_req     ,recv_sta,ierr)
   is = ddc_grid%id+1 ! used as temporary
   if(ddc_grid%id.ne.ide) call mpi_wait(req,recv_sta(:,is),ierr)
   deallocate(recv_sta)

   ! Let us also close the feedback sends
   allocate(recv_sta(mpi_status_size,0:ddc_grid%id-1))
   call mpi_waitall(size(send_conf_req),send_conf_req,recv_sta,ierr)
   deallocate(recv_sta)

   ! And finally get the global values for the counters
   ! Note that gdof already contains the right value
   call mpi_bcast(gdof(1),grid%d,mpi_integer,ide,comm,ierr)
   gdof_dcounter(1:grid%d) = gdof

 end subroutine send_ddc_data

!-----------------------------------------------------------------------
 
 !> Set the global dir and nat numbering
 !!
 !! In this subroutine we don't differentiate vertexes and faces, so
 !! that it's simpler to make a dedicated global communication rather
 !! than including this information in the previous send and receive
 !! operations. The algorithm is composed of three steps: receive data
 !! from lower id subdomains, define the remaining data, send data to
 !! higher id subdomains.
 !<
 subroutine set_g_dir_nat_data(ddc_dofs,dofs,ddc_grid,comm)
  type(t_cgdofs), intent(in) :: dofs
  type(t_ddc_grid), intent(in) :: ddc_grid
  integer, intent(in) :: comm
  type(t_ddc_cgdofs), intent(inout) :: ddc_dofs
 
  integer, parameter :: dim1 = 2
  integer :: is, i, in, j, id, ierr
  !> Counter for the creation of dir and nat dofs. In the creation of
  !! nat dofs, we have to respect the global constraint that element
  !! dofs are numbered after vertex and face dofs. For this reason,
  !! the three counter components count: dir dofs, vertex and face nat
  !! dofs, element nat dofs.
  integer :: dir_nat_cnt(3)
  integer, allocatable :: buf(:,:), req(:), sta(:,:), ioff(:), &
    g_dir_nat_dofs(:)
  character(len=10000) :: msg
  character(len=*), parameter :: &
    this_sub_name = 'set_g_dir_nat_data'

   ! 1) Receive data
   allocate(g_dir_nat_dofs(dofs%ndofs)); g_dir_nat_dofs = -1
   if(ddc_grid%id.ne.0) then
     ! 1.1) MPI communications
     allocate( buf(dim1,sum(ddc_dofs%nndofs_id(ils:ile  ))) )
     allocate( req(                            ils:ile+1  ) )
     allocate( sta( mpi_status_size ,          ils:ile+1  ) )
     req = mpi_request_null
     is = 1
     do i=ils,ile
       if(ddc_dofs%nndofs_id(i).gt.0) &
         call mpi_irecv( buf(1,is) , dim1*ddc_dofs%nndofs_id(i) ,  &
                         mpi_integer, i, 5, comm, req(i), ierr)
       if(i.eq.ddc_grid%id-1) call mpi_irecv( dir_nat_cnt(1) , 3 , &
                         mpi_integer, i, 6, comm, req(ile+1), ierr )
       is = is + ddc_dofs%nndofs_id(i)
     enddo
     call mpi_waitall(size(req),req,sta,ierr)
     ! 1.2) read the received information
     do i=1,size(buf,2)
       in = buf(1,i)
       if(g_dir_nat_dofs(in).eq.-1) then
         g_dir_nat_dofs(in) = buf(2,i)
       else ! consistency check
         if(g_dir_nat_dofs(in).ne.buf(2,i)) then
           write(msg,'(a,i3,a,i8,a,i8)') &
             'Inconsistency on proc ', ddc_grid%id,    &
             '; conflicting global indexes ', buf(2,i), &
             ' and ', g_dir_nat_dofs(in)
           call error(this_sub_name,this_mod_name,msg)
         endif
       endif
     enddo
     deallocate( buf , req , sta )
   else
     dir_nat_cnt = 0
   endif

   ! 2) Define data
   do i=1,size(dofs%dir_dofs,2) ! dir dofs
     in = dofs%dir_dofs(1,i)
     if(g_dir_nat_dofs(in).eq.-1) then
       dir_nat_cnt(1) = dir_nat_cnt(1) + 1
       g_dir_nat_dofs(in) = dir_nat_cnt(1)
     endif
   enddo
   ! The natural dofs are generated in two steps: first the vertex and
   ! face dofs, and then the element dofs.
   do i=1,size(dofs%nat_dofs)-dofs%ndofsd(dofs%d)
     in = dofs%nat_dofs(i)
     if(g_dir_nat_dofs(in).eq.-1) then
       dir_nat_cnt(2) = dir_nat_cnt(2) + 1
       g_dir_nat_dofs(in) = dir_nat_cnt(2)
     endif
   enddo
   do i=size(dofs%nat_dofs)-dofs%ndofsd(dofs%d)+1,size(dofs%nat_dofs)
     in = dofs%nat_dofs(i)
     ! these dofs are never shared
     dir_nat_cnt(3) = dir_nat_cnt(3) + 1
     g_dir_nat_dofs(in) = dir_nat_cnt(3)
   enddo

   ! 3) Send data
   ! 3.1) Prepare the buffer
   allocate( buf(dim1,sum(ddc_dofs%nndofs_id(ihs  :ihe))) )
   allocate( req(                            ihs-1:ihe  ) )
   allocate( sta( mpi_status_size ,          ihs-1:ihe  ) )
   allocate(ioff(ihs:ihe)) ! offset array to fill the buffer
   ioff = 0 ! the size of ioff can be zero
   do i=ihs+1,ihe
     ioff(i) = ioff(i-1) + ddc_dofs%nndofs_id(i-1)
   enddo
   do i=1,ddc_dofs%nndofs
     in = ddc_dofs%ndofs(i)%i
     do j=1,ddc_dofs%ndofs(i)%nd
       id = ddc_dofs%ndofs(i)%id(j)
       if(id.gt.ddc_grid%id) then
         ioff(id) = ioff(id) + 1
         buf(1,ioff(id)) = ddc_dofs%ndofs(i)%in(j)
         buf(2,ioff(id)) = g_dir_nat_dofs(in)
       endif
     enddo
   enddo
   deallocate(ioff)
   ! 3.2) MPI communications
   req = mpi_request_null
   is = 1
   do i=ihs,ihe
     if(ddc_dofs%nndofs_id(i).gt.0) &
       call mpi_isend( buf(1,is) , dim1*ddc_dofs%nndofs_id(i) ,  &
                       mpi_integer, i, 5, comm, req(i), ierr)
     if(i.eq.ddc_grid%id+1) call mpi_isend( dir_nat_cnt(1) , 3 , &
                       mpi_integer, i, 6, comm, req(ihs-1), ierr )
     is = is + ddc_dofs%nndofs_id(i)
   enddo
   call mpi_waitall(size(req),req,sta,ierr)
   deallocate( buf , req , sta )

   ! The last subdomain now broadcasts the final counter values
   call mpi_bcast(dir_nat_cnt(1),3,mpi_integer,ide,comm,ierr)
   ddc_dofs%ngdir_dofs = dir_nat_cnt(1)
   ddc_dofs%ngnat_dofs = sum(dir_nat_cnt(2:3))
   ! fix the global numbers of the element dofs
   do i=size(dofs%nat_dofs)-dofs%ndofsd(dofs%d)+1,size(dofs%nat_dofs)
     in = dofs%nat_dofs(i)
     g_dir_nat_dofs(in) = g_dir_nat_dofs(in) + dir_nat_cnt(2)
   enddo
   allocate(ddc_dofs%gdir_dofs(size(dofs%dir_dofs,2)))
   allocate(ddc_dofs%gnat_dofs(size(dofs%nat_dofs  )))
   ddc_dofs%gdir_dofs = g_dir_nat_dofs(dofs%dir_dofs(1,:))
   ddc_dofs%gnat_dofs = g_dir_nat_dofs(dofs%nat_dofs     )
   deallocate( g_dir_nat_dofs )

 end subroutine set_g_dir_nat_data
 
!-----------------------------------------------------------------------

 subroutine set_ghost_dofs(ddc_dofs,grid,dofs,ddc_grid,comm)
  type(t_grid), intent(in) :: grid
  type(t_cgdofs), intent(in) :: dofs
  type(t_ddc_grid), intent(in) :: ddc_grid
  integer, intent(in) :: comm
  type(t_ddc_cgdofs), intent(inout) :: ddc_dofs

  integer :: ndofs_el, iim1, rim1, i, j, id, is, ie, ierr
  integer, allocatable :: work(:,:), ioff(:), dofs_el(:), &
    iendbuf(:,:), iecvbuf(:,:), added(:,:)
  real(wp), allocatable :: rendbuf(:,:), recvbuf(:,:)

   ndofs_el = size(dofs%dofs,1) ! dofs for each element
   allocate( ddc_dofs%ghost_dofs(ndofs_el,ddc_grid%nns) )
   if(ddc_grid%nns.eq.0) return ! nothing else to do

   ! build the work array
   allocate(work(dofs%ndofs,5))
   work(dofs%nat_dofs     ,1) = 0 ! nat
   work(dofs%nat_dofs     ,2) = ddc_dofs%gnat_dofs
   work(dofs%nat_dofs     ,3) = (/(i, i=1,size(dofs%nat_dofs  ))/)
   work(dofs%dir_dofs(1,:),1) = 1 ! dir
   work(dofs%dir_dofs(1,:),2) = ddc_dofs%gdir_dofs
   work(dofs%dir_dofs(1,:),3) = (/(i, i=1,size(dofs%dir_dofs,2))/)
   work(dofs%dir_dofs(1,:),4) = dofs%dir_dofs(2,:)
   work(dofs%dir_dofs(1,:),5) = dofs%dir_dofs(3,:)

   ! The remaining of this subroutine is similar to
   ! define_ghost_entities in mod_bcs

   ! 1) Communication
   allocate(ioff(0:ddc_grid%nd-1))
   ioff(0) = 0
   do i=1,ddc_grid%nd-1
     ioff(i) = ioff(i-1) + ddc_grid%nns_id(i-1)
   enddo
   iim1 = 1 + 6*ndofs_el
   rim1 = grid%d*ndofs_el
   allocate( iendbuf(iim1,ddc_grid%nns) , iecvbuf(iim1,ddc_grid%nns) ,&
             rendbuf(rim1,ddc_grid%nns) , recvbuf(rim1,ddc_grid%nns) )
   allocate( dofs_el(ndofs_el) )
   do i=1,ddc_grid%nns
     id = ddc_grid%ns(i)%id
     ioff(id) = ioff(id) + 1
     dofs_el = dofs%dofs(:,grid%s(ddc_grid%ns(i)%i)%ie(1)) 

     ! integer buffer
     is = 1; ie = is ! side local position on neighbour
     iendbuf(is,ioff(id)) = ddc_grid%ns(i)%in
     is = ie + 1; ie = is + ndofs_el -1
     call fill_local_idx( iendbuf(is:ie,ioff(id)) ,ddc_dofs,dofs_el,id)
     is = ie + 1; ie = is + ndofs_el -1
     iendbuf(is:ie,ioff(id)) = ddc_dofs%gdofs(dofs_el)%gi
     is = ie + 1; ie = is + ndofs_el - 1
     iendbuf(is:ie,ioff(id)) = work(dofs_el,1) ! dir
     is = ie + 1; ie = is + ndofs_el - 1
     iendbuf(is:ie,ioff(id)) = work(dofs_el,2) ! gi_dir_nat
     is = ie + 1; ie = is + ndofs_el - 1
     iendbuf(is:ie,ioff(id)) = work(dofs_el,4) ! btype
     is = ie + 1; ie = is + ndofs_el - 1
     iendbuf(is:ie,ioff(id)) = work(dofs_el,5) ! breg

     ! real buffer
     do j=1,ndofs_el
       is = (j-1)*grid%d + 1
       ie = is + grid%d - 1
       rendbuf(is:ie,ioff(id)) = dofs%x(:,dofs_el(j))
     enddo
   enddo
   ioff(0) = 0
   do i=1,ddc_grid%nd-1
     ioff(i) = ioff(i-1) + ddc_grid%nns_id(i-1)
   enddo
   call mpi_alltoallv(                                      &
     iendbuf, iim1*ddc_grid%nns_id, iim1*ioff, mpi_integer, &
     iecvbuf, iim1*ddc_grid%nns_id, iim1*ioff, mpi_integer, &
                       comm,ierr )
   call mpi_alltoallv(                                      &
     rendbuf, rim1*ddc_grid%nns_id, rim1*ioff, wp_mpi, &
     recvbuf, rim1*ddc_grid%nns_id, rim1*ioff, wp_mpi, &
                       comm,ierr )

   ! 2) Fill the ghost dofs
   do i=1,ddc_grid%nns
     id = ddc_grid%gs(iecvbuf(1,i))%ni ! local position of the ddc side

     is = 2; ie = is + ndofs_el - 1
     ddc_dofs%ghost_dofs(:,id)%i          = iecvbuf(is:ie,i)
     is = ie + 1; ie = is + ndofs_el - 1
     ddc_dofs%ghost_dofs(:,id)%gi         = iecvbuf(is:ie,i)
     is = ie + 1; ie = is + ndofs_el - 1
     ddc_dofs%ghost_dofs(:,id)%dir        = iecvbuf(is:ie,i).eq.1
     is = ie + 1; ie = is + ndofs_el - 1
     ddc_dofs%ghost_dofs(:,id)%gi_dir_nat = iecvbuf(is:ie,i)
     is = ie + 1; ie = is + ndofs_el - 1
     ddc_dofs%ghost_dofs(:,id)%btype      = iecvbuf(is:ie,i)
     is = ie + 1; ie = is + ndofs_el - 1
     ddc_dofs%ghost_dofs(:,id)%breg       = iecvbuf(is:ie,i)

     do j=1,ndofs_el
       allocate(ddc_dofs%ghost_dofs(j,id)%x(grid%d))
       is = (j-1)*grid%d + 1
       ie = is + grid%d - 1
       ddc_dofs%ghost_dofs(j,id)%x = recvbuf(is:ie,i) 
     enddo
   enddo
   ! Number the ghost dofs which where not know before
   ddc_dofs%n_ghost_dofs = 0
   ddc_dofs%n_ghost_dir_dofs = 0
   ddc_dofs%n_ghost_nat_dofs = 0
   allocate( added(ddc_grid%nns*ndofs_el,2) )
   do i=1,ddc_grid%nns
     do id=1,ndofs_el
       if(ddc_dofs%ghost_dofs(id,i)%i.eq.-1) then
         ddc_dofs%n_ghost_dofs = ddc_dofs%n_ghost_dofs + 1
         ddc_dofs%ghost_dofs(id,i)%i = &
           dofs%ndofs + ddc_dofs%n_ghost_dofs 
         added(ddc_dofs%n_ghost_dofs,:) = (/ id , i /)
         if(ddc_dofs%ghost_dofs(id,i)%dir) then
           ddc_dofs%n_ghost_dir_dofs = ddc_dofs%n_ghost_dir_dofs + 1
           ddc_dofs%ghost_dofs(id,i)%i_dir_nat = &
             size(dofs%dir_dofs,2) + ddc_dofs%n_ghost_dir_dofs 
         else
           ddc_dofs%n_ghost_nat_dofs = ddc_dofs%n_ghost_nat_dofs + 1
           ddc_dofs%ghost_dofs(id,i)%i_dir_nat = &
             size(dofs%nat_dofs  ) + ddc_dofs%n_ghost_nat_dofs 
         endif
       else
         ddc_dofs%ghost_dofs(id,i)%i_dir_nat = &
           work(ddc_dofs%ghost_dofs(id,i)%i,3)
       endif
     enddo
   enddo
 
   deallocate( work , ioff , dofs_el , iendbuf , iecvbuf , &
               rendbuf , recvbuf )

   allocate(ddc_dofs%added_ghost_dofs(ddc_dofs%n_ghost_dofs,2))
   ddc_dofs%added_ghost_dofs = added(1:ddc_dofs%n_ghost_dofs,:)
   deallocate( added )

  contains

   pure subroutine fill_local_idx(idx,ddc_dofs,dofs_el,id)
    type(t_ddc_cgdofs), intent(in) :: ddc_dofs
    integer, intent(in) :: dofs_el(:), id
    integer, intent(out) :: idx(:)

    integer :: i, ni, k

     idx = -1
     do i=1,size(dofs_el)
       if(ddc_dofs%gdofs(dofs_el(i))%ni.gt.0) then
         ! loop on the connected subdomains to see whether the dof is
         ! already defined on neighbour id
         ni = ddc_dofs%gdofs(dofs_el(i))%ni
         check_do: do k=1,ddc_dofs%ndofs(ni)%nd
           if(ddc_dofs%ndofs(ni)%id(k).eq.id) then
             idx(i) = ddc_dofs%ndofs(ni)%in(k)
             exit check_do
           endif
         enddo check_do
       endif
     enddo
   end subroutine fill_local_idx

 end subroutine set_ghost_dofs

!-----------------------------------------------------------------------

 pure subroutine clear_cgdofs( cgdofs )
  type(t_cgdofs), intent(out) :: cgdofs
 end subroutine clear_cgdofs
 
!-----------------------------------------------------------------------

 pure subroutine clear_ddc_cgdofs( ddc_cgdofs )
  type(t_ddc_cgdofs), intent(out) :: ddc_cgdofs
 end subroutine clear_ddc_cgdofs
 
!-----------------------------------------------------------------------

 pure subroutine add_alloc(v,i)
  integer, allocatable, intent(inout) :: v(:)
  integer, intent(in) :: i

  integer :: tmp(size(v)), s

   s = size(v); tmp = v
   deallocate(v); allocate(v(s+1))
   v(1:s) = tmp; v(s+1) = i

 end subroutine add_alloc

!-----------------------------------------------------------------------

 subroutine write_cgdofs(dofs,var_name,fu)
  integer, intent(in) :: fu
  type(t_cgdofs), intent(in) :: dofs
  character(len=*), intent(in) :: var_name
 
  character(len=*), parameter :: &
    this_sub_name = 'write_cgdofs'
 
   write(fu,'(a,a)')    '# name: ',var_name
   write(fu,'(a)')      '# type: struct'
   write(fu,'(a)')      '# length: 9' ! number of fields

   ! field 01 : d
   write(fu,'(a)')      '# name: d'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(dofs%d,'<cell-element>',fu)

   ! field 02 : m
   write(fu,'(a)')      '# name: m'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(dofs%m,'<cell-element>',fu)

   ! field 03 : ndofs
   write(fu,'(a)')      '# name: ndofs'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(dofs%ndofs,'<cell-element>',fu)

   ! field 04 : ndofsd
   write(fu,'(a)')      '# name: ndofsd'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(dofs%ndofsd,'r','<cell-element>',fu)

   ! field 05 : dofs
   write(fu,'(a)')      '# name: dofs'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(dofs%dofs,'<cell-element>',fu)

   ! field 06 : x
   write(fu,'(a)')      '# name: x'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(dofs%x,'<cell-element>',fu)

   ! field 07 : dir_dofs
   write(fu,'(a)')      '# name: dir_dofs'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(dofs%dir_dofs,'<cell-element>',fu)

   ! field 08 : idir_dofs
   write(fu,'(a)')      '# name: idir_dofs'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(dofs%idir_dofs,'<cell-element>',fu)

   ! field 09 : nat_dofs
   write(fu,'(a)')      '# name: nat_dofs'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(dofs%nat_dofs,'r','<cell-element>',fu)

 end subroutine write_cgdofs

!-----------------------------------------------------------------------

 subroutine write_ddc_cgdofs(dofs,var_name,fu)
  integer, intent(in) :: fu
  type(t_ddc_cgdofs), intent(in) :: dofs
  character(len=*), intent(in) :: var_name
 
  integer, allocatable :: tmp(:,:)
  character(len=*), parameter :: &
    this_sub_name = 'write_ddc_cgdofs'
 
   write(fu,'(a,a)')    '# name: ',var_name
   write(fu,'(a)')      '# type: struct'
   write(fu,'(a)')      '# length: 14' ! number of fields

   ! field 01 : ngdofs
   write(fu,'(a)')      '# name: ngdofs'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(dofs%ngdofs,'<cell-element>',fu)

   ! field 02 : ngdofsd
   write(fu,'(a)')      '# name: ngdofsd'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(dofs%ngdofsd,'r','<cell-element>',fu)

   ! field 03 : nndofs
   write(fu,'(a)')      '# name: nndofs'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(dofs%nndofs,'<cell-element>',fu)

   ! field 04 : nndofs_id
   write(fu,'(a)')      '# name: nndofs_id'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(dofs%nndofs_id,'r','<cell-element>',fu)

   ! field 05 : gdofs
   write(fu,'(a)')      '# name: gdofs'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   allocate(tmp(2,size(dofs%gdofs)))
   tmp(1,:) = dofs%gdofs%gi
   tmp(2,:) = dofs%gdofs%ni
   call write_octave(tmp,'<cell-element>',fu)
   deallocate(tmp)

   ! field 06 : ndofs
   write(fu,'(a)')      '# name: ndofs'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(dofs%ndofs,'<cell-element>',fu)

   ! field 07 : ngdir_dofs
   write(fu,'(a)')      '# name: ngdir_dofs'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(dofs%ngdir_dofs,'<cell-element>',fu)

   ! field 08 : ngnat_dofs
   write(fu,'(a)')      '# name: ngnat_dofs'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(dofs%ngnat_dofs,'<cell-element>',fu)

   ! field 09 : gdir_dofs
   write(fu,'(a)')      '# name: gdir_dofs'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(dofs%gdir_dofs,'r','<cell-element>',fu)

   ! field 10 : gnat_dofs
   write(fu,'(a)')      '# name: gnat_dofs'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(dofs%gnat_dofs,'r','<cell-element>',fu)

   ! field 11 : ghost_dofs
   write(fu,'(a)')      '# name: ghost_dofs'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(dofs%ghost_dofs,'<cell-element>',fu)

   ! field 12 : n_ghost_dofs
   write(fu,'(a)')      '# name: n_ghost_dofs'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(dofs%n_ghost_dofs,'<cell-element>',fu)

   ! field 13 : n_ghost_dir_dofs
   write(fu,'(a)')      '# name: n_ghost_dir_dofs'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(dofs%n_ghost_dir_dofs,'<cell-element>',fu)

   ! field 14 : n_ghost_nat_dofs
   write(fu,'(a)')      '# name: n_ghost_nat_dofs'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(dofs%n_ghost_nat_dofs,'<cell-element>',fu)

 end subroutine write_ddc_cgdofs

!-----------------------------------------------------------------------

 subroutine write_ghost_dof_struct(g,var_name,fu)
  integer, intent(in) :: fu
  type(t_ghost_dof), intent(in) :: g(:,:)
  character(len=*), intent(in) :: var_name
 
  integer :: i, j
  character(len=*), parameter :: &
    this_sub_name = 'write_ghost_dof_struct'
 
   write(fu,'(a,a)')    '# name: ',var_name
   write(fu,'(a)')      '# type: struct'
   write(fu,'(a)')      '# length: 8' ! number of fields

   ! field 01 : i
   write(fu,'(a)')      '# name: i'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a,i7)')   '# rows: ',size(g,1)
   write(fu,'(a,i7)')   '# columns: ',size(g,2)
   do j=1,size(g,2)
     do i=1,size(g,1)
       call write_octave(g(i,j)%i,'<cell-element>',fu)
     enddo
   enddo

   ! field 02 : i_dir_nat
   write(fu,'(a)')      '# name: i_dir_nat'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a,i7)')   '# rows: ',size(g,1)
   write(fu,'(a,i7)')   '# columns: ',size(g,2)
   do j=1,size(g,2)
     do i=1,size(g,1)
       call write_octave(g(i,j)%i_dir_nat,'<cell-element>',fu)
     enddo
   enddo

   ! field 03 : gi
   write(fu,'(a)')      '# name: gi'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a,i7)')   '# rows: ',size(g,1)
   write(fu,'(a,i7)')   '# columns: ',size(g,2)
   do j=1,size(g,2)
     do i=1,size(g,1)
       call write_octave(g(i,j)%gi,'<cell-element>',fu)
     enddo
   enddo

   ! field 04 : dir
   write(fu,'(a)')      '# name: dir'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a,i7)')   '# rows: ',size(g,1)
   write(fu,'(a,i7)')   '# columns: ',size(g,2)
   do j=1,size(g,2)
     do i=1,size(g,1)
       call write_octave(g(i,j)%dir,'<cell-element>',fu)
     enddo
   enddo

   ! field 05 : gi_dir_nat
   write(fu,'(a)')      '# name: gi_dir_nat'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a,i7)')   '# rows: ',size(g,1)
   write(fu,'(a,i7)')   '# columns: ',size(g,2)
   do j=1,size(g,2)
     do i=1,size(g,1)
       call write_octave(g(i,j)%gi_dir_nat,'<cell-element>',fu)
     enddo
   enddo

   ! field 06 : x
   write(fu,'(a)')      '# name: x'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a,i7)')   '# rows: ',size(g,1)
   write(fu,'(a,i7)')   '# columns: ',size(g,2)
   do j=1,size(g,2)
     do i=1,size(g,1)
       call write_octave(g(i,j)%x,'c','<cell-element>',fu)
     enddo
   enddo

   ! field 07 : btype
   write(fu,'(a)')      '# name: btype'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a,i7)')   '# rows: ',size(g,1)
   write(fu,'(a,i7)')   '# columns: ',size(g,2)
   do j=1,size(g,2)
     do i=1,size(g,1)
       call write_octave(g(i,j)%btype,'<cell-element>',fu)
     enddo
   enddo

   ! field 08 : breg
   write(fu,'(a)')      '# name: breg'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a,i7)')   '# rows: ',size(g,1)
   write(fu,'(a,i7)')   '# columns: ',size(g,2)
   do j=1,size(g,2)
     do i=1,size(g,1)
       call write_octave(g(i,j)%breg,'<cell-element>',fu)
     enddo
   enddo

 end subroutine write_ghost_dof_struct
 
!-----------------------------------------------------------------------

end module mod_cgdofs

